clear;clc;close all;filebrowser;workspace;warning('off','all');
%% Set flags
massFlag = 1;
%%-----MCP Parameters--------------
TMCPacquire = 2*10^5; %set in the MCP acquisition software, in [ns]
dtMCPres = 0.025117; %TDC resolution [ns]
MCPRadius = 42;     %MCP detector radius [mm]
%%-----Parameters for data analysis--------
Xc = -10.250*0;      %[mm] center for plotting in detector frame
Yc = 8.75*0;    %[mm] center for plotting in detector frame
Angle = -30*0;   %[deg] rotation angle for plotting 
tbinsize = 5; %20 [ns] ignore this if you are histcounts mass
Nbin = floor(TMCPacquire/tbinsize);
tEdge = (0:1:Nbin)*tbinsize; %[ns] for TOF binning
tList = (0.5:1:(Nbin-0.5))*tbinsize/1000; %[us]
% f(us) = 0.0214+ 2.481*sqrt(mass(amu));  %This formula from Yu's calibration
% massList = (tList/1000./15.69).^2*40-0.26;  %[amu] mass
if massFlag
    dmass = 0.25;
    massEdge = (20-dmass/2):dmass:(300+dmass/2);    %[amu]
    massList = 20:dmass:300;        %[amu]
%     tList = 21.4 + 2481*sqrt(massList);  %[ns]
end
dXY = 0.5;        %[mm] spatial resolution
Xmin = -45;     %[mm] plot range
Ymin = Xmin;    %[mm]
Xmax = 45;      %[mm]
Ymax = Xmax;    %[mm]
Xedge = Xmin:dXY:Xmax;
Yedge = Ymin:dXY:Ymax;
Xplot = (Xmin+dXY/2):dXY:(Xmax-dXY/2);
Yplot = (Ymin+dXY/2):dXY:(Ymax-dXY/2);
LabelFontSize = 12;
dmTOF = 6;      %% mass range for removing dominant peak in mass spectrum
mlist1 = [40 87 127];  
species1 = {'K^+' 'Rb^+' 'KRb^+'};
dmPlot = 6;     %% mass range for plotting [amu]
mlistPlot = [127 87 40 80 174 214 167 254];%for spatial plot
speciesPlot = {'KRb^+' 'Rb^+' 'K^+' 'K_2^+' 'Rb_2^+' 'KRb_2^+' 'K_2Rb^+' 'K_2Rb_2^+'};%for spatial plot
colorbg = zeros(1,3);       %Black
% colorbg = ones(1,3);        %white
%%---------Define folderpath-----------
% folderpath = 'C:\Users\KRb\lmf2txt';
% folderpath = '/Users/MGHu/Documents/KRb reaction/Data/Figure 3 data';
%% ------load .dat files--------------------------
% %% solid beam data 1
% filenameList={'2018_9_27_6'};
% NexpCycle = [floor(2213/1)]; % Not MCP trigger
% t0 = 15;     %[ms]
%% Solid beam data 2
% filenameList = {'2018_11_28_3_1' '2018_11_28_3_2' '2018_11_28_4' '2018_12_3_1'};
% NexpCycle = [1649 floor((15234-1649-(11176-9106))/10) floor(12416/10) floor(42216/30)]; % Not MCP trigger
% t0 = 15;     %[ms]

% %% Hollow-bottle beam data 1---
% filenameList = {'2018_12_7_1'};
% NexpCycle = [floor(107516/30)]; % Not MCP trigger

% %% Hollow-bottle beam data 2---
% filenameList = {'2018_12_10_4'};
% NexpCycle = [floor(20997/30)]; % Not MCP trigger
% t0 = 50;     %[ms]

% %% cross Hollow-bottle beam data 1(see readme.txt for detail)---
% filenameList = {'2019_01_03_2'};
% NexpCycle = [floor(124409/30)]; % Not MCP trigger
% t0 = 30;     %[ms]  time delay between STIRAP and 1st UV pulse

% %% cross Hollow-bottle beam background data of Rb atoms (see readme.txt for detail)---
% filenameList = {'2019_01_07_1'};
% NexpCycle = [floor(124394/100)]; % Not MCP trigger
% t0 = 30;     %[ms]  time delay between STIRAP and 1st UV pulse
 
% %% cross Hollow-bottle beam background data (see readme.txt for detail)---
% filenameList = {'2019_01_03_1'};
% NexpCycle = [floor(4995/100)]; % Not MCP trigger
% t0 = 15;     %[ms]  time delay between STIRAP and 1st UV pulse

% %% cross Hollow-bottle beam data 296nm(see readme.txt for detail)---
% filenameList = {'2019_01_12_7'};
% NexpCycle = [floor(83154/20)]; % Not MCP trigger
% t0 = 40;     %[ms]  time delay between STIRAP and 1st UV pulse

% %% cross Hollow-bottle beam data 300nm(see readme.txt for detail)---
% filenameList = {'2019_01_21_1'};
% NexpCycle = [floor(215840/20)]; % Not MCP trigger
% t0 = 40;     %[ms]  time delay between STIRAP and 1st UV pulse
% 

% %% cross Hollow-bottle beam data 285nm(see readme.txt for detail)---
filenameList = {'2019_01_29_1'};
NexpCycle = [floor(215840/20)]; % Not MCP trigger
t0 = 40;     %[ms]  time delay between STIRAP and 1st UV pulse

% %% cross Hollow-bottle beam data 285nm, longer cleanup pulse(see readme.txt for detail)---
% filenameList = {'2019_02_04_1'};
% NexpCycle = [floor(198022/20)]; % Not MCP trigger
% t0 = 40;     %[ms]  time delay between STIRAP and 1st UV pulse

%% figures
h1 = figure('units','normalized','WindowStyle','docked','outerposition',[0 0 1 1]);
h2 = figure('units','normalized','WindowStyle','docked','outerposition',[0 0 1 1]);
h3 = figure('units','normalized','WindowStyle','docked','outerposition',[0 0 1 1]);
h4 = figure('units','normalized','WindowStyle','docked','outerposition',[0 0 1 1]);
h5 = figure('units','normalized','WindowStyle','docked','outerposition',[0 0 1 1]);
h6 = figure('units','normalized','WindowStyle','docked','outerposition',[0 0 1 1]);
h7 = figure('units','normalized','WindowStyle','docked','outerposition',[0 0 1 1]);

Nfile = length(filenameList);
realHitList = cell(1, Nfile);
NUVpulse = zeros(1,Nfile);
dataTOFList = cell(1, Nfile);
dataMassList = cell(1, Nfile);
dataXYAllList = cell(1, Nfile);
dataXYPlotList = cell(1, Nfile);
for i = 1:Nfile
    filename = filenameList{i};
    if (~exist([filename,'.dat'],'file'))
        error(['Please use mainMCPtxtAnalysis.m to convert /lmf2txt/',filename,...
            '.lmf.txt to ' filename, '.dat']);
    else
        realHitList{i} = dlmread([filename,'.dat'],'\t'); 
        %1st column ---- TOF [ns]
        %2nd column ---- X position [mm]
        %3rd column ---- Y position [mm]
        %4th column ---- MCP trigger #
        %5th column ---- Exp cycle #
    end
    NUVpulse(i) = max(realHitList{i}(:,4));
    dataTOF = zeros(NUVpulse(i),length(tList)); 
    dataMass = zeros(NUVpulse(i), length(massList));
    xtAll = realHitList{i}(:,1);
    XAll = realHitList{i}(:,2) - Xc;
    YAll = realHitList{i}(:,3) - Yc;
    dataXYAll = cell(1, NUVpulse(i));
    dataXYi = cell(1, NUVpulse(i));
    for j = 1:NUVpulse(i)
        xt = xtAll(realHitList{i}(:,4) == j); %[ns]
        dataTOF(j,:)= histcounts(xt, tEdge);
        xmass = ((xt-21.4)/2481).^2;  %[amu]
        dataMass(j,:) = histcounts(xmass, massEdge);
        %--All species at UVpulse j----
        dataX0 = XAll(realHitList{i}(:,4) == j);
        dataY0 = YAll(realHitList{i}(:,4) == j);
        dataX = dataX0*cos(-Angle*pi/180)-dataY0*sin(-Angle*pi/180);
        dataY = dataX0*sin(-Angle*pi/180)+dataY0*cos(-Angle*pi/180);
        dataXYAll{j} = histcounts2(dataY,dataX,Yedge,Xedge);
        dataXYj = zeros(length(Yplot),length(Xplot),length(mlistPlot));
        %--selected species at UVpulse j--------
        for k = 1:length(mlistPlot)
            dataRbX = dataX((xmass>=(mlistPlot(k)-dmPlot/2))&(xmass<=(mlistPlot(k)+dmPlot/2)));
            dataRbY = dataY((xmass>=(mlistPlot(k)-dmPlot/2))&(xmass<=(mlistPlot(k)+dmPlot/2)));
            dataXYj(:,:,k) = histcounts2(dataRbY,dataRbX,Yedge,Xedge);
        end
        dataXYi{j} = dataXYj;
    end
    dataTOFList{i} = dataTOF;       %[ns]
    dataXYAllList{i} = dataXYAll;      %[mm]
    dataXYPlotList{i} = dataXYi;
    dataMassList{i} = dataMass;     %[amu]    
end
disp('UV pulse per cycle are:');
disp(NUVpulse);
disp('Experiment cycle is');
disp(num2str(NexpCycle));

%%-----Plot combined density-dependent mass spectrum-----------------------
if 1
    figure(h2);
    UVpulseNoList = [1 2 3 4 5 6];
    for k = 1:length(UVpulseNoList)
        UVpulseNo = UVpulseNoList(k);
        subplot(2,3,k)
        if length(NexpCycle)>1
            if UVpulseNo < 2
                dataMassSum = dataMassList{1}(UVpulseNo,:)+dataMassList{2}(UVpulseNo,:)+dataMassList{3}(UVpulseNo,:)+dataMassList{4}(UVpulseNo,:);
            elseif UVpulseNo <11
                dataMassSum = (dataMassList{2}(UVpulseNo,:)+dataMassList{3}(UVpulseNo,:)+dataMassList{4}(UVpulseNo,:))*sum(NexpCycle)/sum(NexpCycle(2:length(NexpCycle)));
            else
                dataMassSum = dataMassList{4}(UVpulseNo,:)*sum(NexpCycle)/sum(NexpCycle(4:length(NexpCycle)));
            end
        else
            dataMassSum = dataMassList{1}(UVpulseNo,:);
        end
        mTOF = massList;
        countsTOF1 = dataMassSum;
        for i = 1:length(mlist1)
            countsTOF1(mTOF>=(mlist1(i)-dmTOF/2)& mTOF<=(mlist1(i)+dmTOF/2)) = 0;
        end
        plot(massList, countsTOF1,'Color','black');
        hold on
        xlim([20,300]);
        ylim([0,40]);
        mlist = [40 87 127];
        for i = 1:length(mlist)
            xdata = mTOF(mTOF>=(mlist(i)-dmTOF/2)&mTOF<=(mlist(i)+dmTOF/2));
            ydata = dataMassSum(mTOF>=(mlist(i)-dmTOF/2)&mTOF<=(mlist(i)+dmTOF/2));
            plot(xdata, ydata,'-','LineWidth',1, 'Color', [0 0 0]+0.8);
            hold on;
        end
        %     dmTOF = 11;
        mlist = [80 174 214 254 167];
        species = {'K_2^+' 'Rb_2^+' 'KRb_2^+' 'K_2Rb_2^+' 'K_2Rb'};
        colorlist = {'red' 'cyan' 'magenta' 'blue' 'green'};
        for i = 1:length(mlist)
            xdata = mTOF(mTOF>=(mlist(i)-dmTOF/2)&mTOF<=(mlist(i)+dmTOF/2));
            ydata = dataMassSum(mTOF>=(mlist(i)-dmTOF/2)&mTOF<=(mlist(i)+dmTOF/2));
            plot(xdata, ydata,'-','LineWidth', 1.5, 'Color', colorlist{i});
            if strcmp(species{i},'K_2Rb')
                text(mlist(i)-22, 5+max(ydata),species{i}, 'Color',colorlist{i});
            else
                text(mlist(i), 5+max(ydata),species{i}, 'Color',colorlist{i});
            end
            hold on;
        end
        hold off
        title(['UV pulse time after STIRAP: ', num2str((UVpulseNoList(k)-1)*100+t0), ' ms']);
        xlabel('Mass (amu)');
        ylabel('Ion Signal');
    end
end
%%-----Plot combined density-dependent XY distribution of selected species----
if 0
    for i = 1:length(mlistPlot)
        figure('units','normalized','WindowStyle','docked','outerposition',[0 0 1 1]);
        maxCounts = 0;
        for k = 1:length(UVpulseNoList)
            UVpulseNo = UVpulseNoList(k);
            ax = subplot(2,3,k);
            if length(NexpCycle) > 1
                if UVpulseNo < 2
                    dataXYPlot = dataXYPlotList{1}{UVpulseNo}(:,:,i)+dataXYPlotList{2}{UVpulseNo}(:,:,i)+dataXYPlotList{3}{UVpulseNo}(:,:,i)+dataXYPlotList{4}{UVpulseNo}(:,:,i);
                elseif UVpulseNo < 11
                    dataXYPlot = (dataXYPlotList{2}{UVpulseNo}(:,:,i)+dataXYPlotList{3}{UVpulseNo}(:,:,i)+dataXYPlotList{4}{UVpulseNo}(:,:,i))*sum(NexpCycle)/sum(NexpCycle(2:length(NexpCycle)));
                else
                    dataXYPlot = dataXYPlotList{4}{UVpulseNo}(:,:,i)*sum(NexpCycle)/sum(NexpCycle(4:length(NexpCycle)));
                end
            else
                dataXYPlot = dataXYPlotList{1}{UVpulseNo}(:,:,i);
            end
            %     dataXYPlot = dataXYRbList{2}{1};
            imagesc(Xplot, Yplot, imrotate(dataXYPlot, Angle));
            % Make values 0-5 black:
            cmap = jet(1+2*ceil(max(max(dataXYPlot)))); %jet(Number of color)
            if k == 1
                maxCounts = ceil(max(max(dataXYPlot)));
                cmap = jet(1+2*maxCounts);
            end
            if ceil(max(max(dataXYPlot)))>maxCounts
                maxCounts = ceil(max(max(dataXYPlot)));
                cmap = jet(1+2*maxCounts);
            end
            cmap(1:1,:) = colorbg;
            colormap(ax, cmap);
            colorbar
            xlabel('x (mm)');
            ylabel('y (mm)');
            title([speciesPlot{i},': ', num2str((UVpulseNoList(k)-1)*100+t0), ' ms, m=',num2str(mlistPlot(i)), '\pm', num2str(dmPlot/2), ' amu']);
            set(gca,'YDir','normal');
            xtick = get(gca, 'XTick');
            set(gca, 'FontSize', LabelFontSize);
            ytick = get(gca, 'YTick');
            set(gca, 'FontSize', LabelFontSize);
        end
    end
end
%%-----Plot All XY distribution of selected species--------------
if 1
    figure(h3);
    %     dataXYPlot = dataXYRbList{2}{1};
    for i = 1:length(mlistPlot)
        dataXYPlot = zeros(length(Yplot),length(Xplot));
        for j = 1:Nfile
            for k = 1:10%NUVpulse(j)
                dataXYPlot = dataXYPlot+dataXYPlotList{j}{k}(:,:,i);
            end
        end
        subplot(3,3,i)
        maxCounts = ceil(max(max(dataXYPlot)));
        dataXYPlotCut = dataXYPlot;
        if i==2
            dataXYPlotCut(dataXYPlotCut > maxCounts/160) = maxCounts/160;
        end
%         if i==3
%             dataXYPlotCut(dataXYPlotCut > maxCounts/30) = maxCounts/30;
%         end
        if maxCounts < 2
            dataXYPlotCut(1,1) = 2;
        end
        imagesc(Xplot, Yplot, dataXYPlotCut);
        if i == 1    
            cmap = jet(maxCounts);
        end
        cmap(1:1,:) = colorbg;    % Make values 0 black:
        colormap(cmap);
        colorbar
        xlabel('x (mm)');
        ylabel('y (mm)');
        title(['All ',speciesPlot{i},': m=',num2str(mlistPlot(i)), '\pm', num2str(dmPlot/2), ' amu']);
        set(gca,'YDir','normal');
        xtick = get(gca, 'XTick');
        set(gca, 'FontSize', LabelFontSize);
        ytick = get(gca, 'YTick');
        set(gca, 'FontSize', LabelFontSize);
        hold on
        t = linspace(0,2*pi);
        plot(MCPRadius.*cos(t),MCPRadius.*sin(t),'w--');
    end
end
%%----Plot total Mass spectrometer-----------------
if 1
    figure(h1);
    dataMassSum = zeros(size(dataMassList{1}(1,:)));
%     for i = 1:1
%         dataMassSum = dataMassSum+dataMassList{1}(i,:);
%     end
    for i=1:Nfile
        for j=1:3%NUVpulse(i)
            dataMassSum = dataMassSum+dataMassList{i}(j,:);
        end
    end
    mTOF = massList;
    countsTOF1 = dataMassSum;
    subplot(2,1,1)
    plot(massList, countsTOF1,'Color','black');
    xlim([20,300]);
    xlabel('Mass (amu)');
    ylabel('Ion Signal');
%     title(['Mass spectrum of all shots']);
    xt = get(gca, 'XTick');
    set(gca, 'FontSize', LabelFontSize);
    xt = get(gca, 'YTick');
    set(gca, 'FontSize', LabelFontSize);
    subplot(2,1,2)
    for i = 1:length(mlist1)
        countsTOF1(mTOF>=(mlist1(i)-dmTOF/2)& mTOF<=(mlist1(i)+dmTOF/2)) = 0;
    end
    plot(massList, countsTOF1,'Color','black');
    hold on
    xlim([20,300]);
    ylim([0,60]);
    mlist = [40 87 127];
    for i = 1:length(mlist)
        xdata = mTOF(mTOF>=(mlist(i)-dmTOF/2)&mTOF<=(mlist(i)+dmTOF/2));
        ydata = dataMassSum(mTOF>=(mlist(i)-dmTOF/2)&mTOF<=(mlist(i)+dmTOF/2));
        plot(xdata, ydata, '-', 'LineWidth', 0.5, 'Color', [0 0 0]+0.8);
        hold on;
    end
    %     dmTOF = 11;
    mlist = [80 174 214 254 167];
    species = {'K_2^+' 'Rb_2^+' 'KRb_2^+' 'K_2Rb_2^+' 'K_2Rb'};
    colorlist = {'red' 'cyan' 'magenta' 'blue' 'green'};
    for i = 1:length(mlist)
        xdata = mTOF(mTOF>=(mlist(i)-dmTOF/2)&mTOF<=(mlist(i)+dmTOF/2));
        ydata = dataMassSum(mTOF>=(mlist(i)-dmTOF/2)&mTOF<=(mlist(i)+dmTOF/2));
        plot(xdata, ydata, '-', 'LineWidth', 0.75, 'Color', colorlist{i});
%         if strcmp(species{i},'K_2Rb^+')
%             text(mlist(i)-22, 5+max(ydata),species{i}, 'Color',colorlist{i});
%         else
%             text(mlist(i), 5+max(ydata),species{i}, 'Color',colorlist{i});
%         end
        hold on;
        disp(['Numb of ', species{i},'=',num2str(sum(ydata))]);
    end
    hold off
    %     title(['UV pulse time after STIRAP: ', num2str((UVpulseNoList(k)-1)*100+15), ' ms']);
    xlabel('Mass (amu)');
    ylabel('Ion Signal');
    xt = get(gca, 'XTick');
    set(gca, 'FontSize', LabelFontSize);
    xt = get(gca, 'YTick');
    set(gca, 'FontSize', LabelFontSize);
end

%%-----plot selected species--------
figure(h4);
i = 4;
dataXYPlot = zeros(length(Yplot),length(Xplot));
for j = 1:Nfile
    for k = 1:3%NUVpulse(j)
        dataXYPlot = dataXYPlot+dataXYPlotList{j}{k}(:,:,i);
    end
end
maxCounts = ceil(max(max(dataXYPlot)));
dataXYPlotCut = dataXYPlot;
if i==2
    dataXYPlotCut(dataXYPlotCut > maxCounts/200) = maxCounts/200;
end
if i==3
    dataXYPlotCut(dataXYPlotCut > maxCounts/20) = maxCounts/20;
end
if i==4
    dataXYPlotCut(1,1) = 6;
end
if maxCounts < 2
    dataXYPlotCut(1,1) = 2;
end
imagesc(Xplot, Yplot, dataXYPlotCut);
if i == 4
    xplc = -8.867;
    yplc = 7.36;
elseif i == 5
        xplc = -10.25;
        yplc = 8.25;
end

dx = 10;
xlim([xplc-dx/2,xplc+dx/2]);
ylim([yplc-dx/2,yplc+dx/2]);
cmap(1:1,:) = colorbg;    % Make values 0 black:
colormap(cmap);
colorbar
xlabel('x (mm)');
ylabel('y (mm)');
% title([speciesPlot{i},': m=',num2str(mlistPlot(i)), '\pm', num2str(dmPlot/2), ' amu']);
set(gca,'YDir','normal');
xtick = get(gca, 'XTick');
set(gca, 'FontSize', LabelFontSize);
ytick = get(gca, 'YTick');
set(gca, 'FontSize', LabelFontSize);
hold on
t = linspace(0,2*pi);
KE = 1.23984;       % [meV] 10 cm^-1
VR = 1000;          % [V] repeller voltage
AVMI = 46.7;        %[mm/sqrt(meV/V)]
R10cm = AVMI*sqrt(KE/VR);
plot(R10cm.*cos(t)+xplc, R10cm.*sin(t)+yplc, 'y--','LineWidth', 1.5);
hold off


%%----Plot total Mass spectrometer with and without KRb -----------------
if 1
    figure(h5);
    dataMassSum = zeros(size(dataMassList{1}(1,:)));
    
    subplot(2,1,1)
    for i=1:Nfile
        for j=1:3
            dataMassSum = dataMassSum+dataMassList{i}(j,:);
        end
    end
    mTOF = massList;
    countsTOF1 = dataMassSum;

    for i = 1:length(mlist1)
        countsTOF1(mTOF>=(mlist1(i)-dmTOF/2)& mTOF<=(mlist1(i)+dmTOF/2)) = 0;
    end
    plot(massList, countsTOF1,'Color','black');
    hold on
    xlim([20,270]);
    ylim([0, 50]);
    mlist = [40 87 127];
    for i = 1:length(mlist)
        xdata = mTOF(mTOF>=(mlist(i)-dmTOF/2)&mTOF<=(mlist(i)+dmTOF/2));
        ydata = dataMassSum(mTOF>=(mlist(i)-dmTOF/2)&mTOF<=(mlist(i)+dmTOF/2));
        plot(xdata, ydata,'-','LineWidth',1, 'Color', [0 0 0]+0.8);
        hold on;
    end
    %     dmTOF = 11;
    mlist = [80 174 214 167];
    species = {'K_2^+' 'Rb_2^+' 'KRb_2^+' 'K_2Rb'};
    colorlist = {'red' 'cyan' 'magenta' 'green'};
    for i = 1:length(mlist)
        xdata = mTOF(mTOF>=(mlist(i)-dmTOF/2)&mTOF<=(mlist(i)+dmTOF/2));
        ydata = dataMassSum(mTOF>=(mlist(i)-dmTOF/2)&mTOF<=(mlist(i)+dmTOF/2));
        plot(xdata, ydata,'-','LineWidth', 1.0, 'Color', colorlist{i});
%         if strcmp(species{i},'K_2Rb')
%             text(mlist(i)-22, 5+max(ydata),species{i}, 'Color',colorlist{i});
%         else
%             text(mlist(i), 5+max(ydata),species{i}, 'Color',colorlist{i});
%         end
        hold on;
        disp(['Numb of ', species{i},'=',num2str(sum(ydata))]);
    end
    hold off
    %     title(['UV pulse time after STIRAP: ', num2str((UVpulseNoList(k)-1)*100+15), ' ms']);
    xlabel('m/Z');
    ylabel('Ion (count)');
    xt = get(gca, 'XTick');
    set(gca, 'FontSize', LabelFontSize);
    xt = get(gca, 'YTick');
    set(gca, 'FontSize', LabelFontSize);
    
    subplot(2,1,2)
    dataMassSum = zeros(size(dataMassList{1}(1,:)));
    for i=1:Nfile
        for j=18:20%NUVpulse(i)
            dataMassSum = dataMassSum+dataMassList{i}(j,:);
        end
    end

    countsTOF1 = dataMassSum;  

    for i = 1:length(mlist1)
        countsTOF1(mTOF>=(mlist1(i)-dmTOF/2)& mTOF<=(mlist1(i)+dmTOF/2)) = 0;
    end
    plot(massList, countsTOF1,'Color','black');
    hold on
    xlim([20, 270]);
    ylim([0, 50]);
    mlist = [40 87 127];
    for i = 1:length(mlist)
        xdata = mTOF(mTOF>=(mlist(i)-dmTOF/2)&mTOF<=(mlist(i)+dmTOF/2));
        ydata = dataMassSum(mTOF>=(mlist(i)-dmTOF/2)&mTOF<=(mlist(i)+dmTOF/2));
        plot(xdata, ydata,'-','LineWidth',1, 'Color', [0 0 0]+0.8);
        hold on;
    end
    %     dmTOF = 11;
    mlist = [80 174 214 167];
    species = {'K_2^+' 'Rb_2^+' 'KRb_2^+' 'K_2Rb'};
    colorlist = {'red' 'cyan' 'magenta' 'green'};
    for i = 1:length(mlist)
        xdata = mTOF(mTOF>=(mlist(i)-dmTOF/2)&mTOF<=(mlist(i)+dmTOF/2));
        ydata = dataMassSum(mTOF>=(mlist(i)-dmTOF/2)&mTOF<=(mlist(i)+dmTOF/2));
        plot(xdata, ydata,'-','LineWidth', 1.0, 'Color', colorlist{i});
%         if strcmp(species{i},'K_2Rb')
%             text(mlist(i)-22, 5+max(ydata),species{i}, 'Color',colorlist{i});
%         else
%             text(mlist(i), 5+max(ydata),species{i}, 'Color',colorlist{i});
%         end
        hold on;
        disp(['Numb of ', species{i},'=',num2str(sum(ydata))]);
    end
    hold off
    %     title(['UV pulse time after STIRAP: ', num2str((UVpulseNoList(k)-1)*100+15), ' ms']);
    xlabel('m/Z');
    ylabel('Ion (count)');
    xt = get(gca, 'XTick');
    set(gca, 'FontSize', LabelFontSize);
    xt = get(gca, 'YTick');
    set(gca, 'FontSize', LabelFontSize);
end


%%----Plot total TOF spectrometer-----------------
if 1
    figure(h6);
    dataTOFSum = zeros(size(dataTOFList{1}(1,:)));

    for i=1:Nfile
        for j=1:3%NUVpulse(i)
            dataTOFSum = dataTOFSum+dataTOFList{i}(j,:);
        end
    end
    mTOF = tList;
    countsTOF1 = dataTOFSum;
    subplot(2,1,1)
    plot(tList, countsTOF1,'Color','black');
    xlim([11,35]);
    xlabel('TOF (us)');
    ylabel('Ion Count');
    title(['TOF spectrum of ions']);
    xt = get(gca, 'XTick');
    set(gca, 'FontSize', LabelFontSize);
    xt = get(gca, 'YTick');
    set(gca, 'FontSize', LabelFontSize);
    subplot(2,1,2)
    dataMassSum = zeros(size(dataMassList{1}(1,:)));
    for i=1:Nfile
        for j=1:3%NUVpulse(i)
            dataMassSum = dataMassSum+dataMassList{i}(j,:);
        end
    end
    countsTOF1 = dataMassSum;
    plot(massList, countsTOF1,'Color','black');
    xlim([20,300]);
    xlabel('Mass (amu)');
    ylabel('Ion Signal');
%     title(['Mass spectrum of ions']);
    xt = get(gca, 'XTick');
    set(gca, 'FontSize', LabelFontSize);
    xt = get(gca, 'YTick');
    set(gca, 'FontSize', LabelFontSize);
end

%%-----Plot All XY distribution of selected species--------------
if 1
    figure(h7);
    dataTOFSum = zeros(size(dataTOFList{1}(1,:)));
    for k= 1:Nfile
        for j=1:3 %NUVpulse(i)
            dataTOFSum = dataTOFSum+dataTOFList{k}(j,:);
        end
    end
j=1;
    for i = [1 5 6]
        dTOF =0.5; % [us]
        xTOF = (sqrt(mlistPlot(i))*2481+21.4)/1000;   %[us]
%             xmass = ((xt-21.4)/2481).^2;  %[amu]
        xdata = tList(tList>=(xTOF-dTOF/2)&tList<=(xTOF+dTOF/2));%[us]
        ydata = dataTOFSum(tList>=(xTOF-dTOF/2)&tList<=(xTOF+dTOF/2));
        if i == 5
            ymaxRb2 = max(ydata);
            xdataRb2 = xTOF;
        end
        if i == 1
            ymaxKRb = max(ydata);
            xdataKRb = xTOF;
        end
        if i == 6
            ymaxKRb2 = max(ydata);
            xdataKRb2 = xTOF;
        end
        subplot(2,2,j)
        plot(xdata, ydata,'o-','LineWidth', 1.0);
        if i==1
            ax = gca;
            ax.YAxis.Exponent = 3;
        end
        
%         title([speciesPlot{i}]);
        xlim([xTOF-dTOF/2, xTOF+dTOF/2]);
        xlabel('TOF (\mus)');
        ylabel('Ion (count)');
        xticks([xTOF-dTOF/2, xTOF, xTOF+dTOF/2])
        xticklabels([xTOF-dTOF/2, xTOF, xTOF+dTOF/2])
        xt = get(gca, 'XTick');
        set(gca, 'FontSize', LabelFontSize);
        xt = get(gca, 'YTick');
        set(gca, 'FontSize', LabelFontSize);
        j = j+1;
    end
%     xdata_Ring = 1:dTOF*1000; %[ns]
%     ydata_Ring = zeros(size(xdata_Ring));
%     ydata_Ring1 = dlmread(['RingBeamConvolve2.dat'],'\t');
%     ydata_Ring1 = ydata_Ring1/max(ydata_Ring1)*9;%*ymaxRb2;
%     xdata_Ring = xdata_Ring + (xdataRb2-dTOF/2)*1000; %[ns]
%     ydata_Ring(xdata_Ring>32685&xdata_Ring<=32685+length(ydata_Ring1)) = ydata_Ring1;
% %     xdata_Ring = ((1:length(ydata_Ring))-length(ydata_Ring)/2)/1000+32.78;  %[us]
%     subplot(2,2,2)
%     hold on
%     plot(xdata_Ring/1000, ydata_Ring,'-r','LineWidth', 2.0);
%     ylim([0, 12]);
%     hold off
%     
%     xdata_Solid = 1:dTOF*1000; %[ns]
%     ydata_Solid = zeros(size(xdata_Solid));
%     ydata_Solid1 = dlmread(['SolidBeam.dat'],'\t');
%     ydata_Solid1 = ydata_Solid1/max(ydata_Solid1)*ymaxKRb;
%     xdata_Solid = xdata_Solid + (xdataKRb-dTOF/2)*1000; %[ns]
%     ydata_Solid(xdata_Solid>27915&xdata_Solid<=27915+length(ydata_Solid1)) = ydata_Solid1;
%     subplot(2,2,1)
%     hold on
%     plot(xdata_Solid/1000, ydata_Solid,'-r','LineWidth', 2.0);
%     hold off
%     subplot(2,2,3)
%     hold on
%     dTOF1 =0.01; % [ns] set by KRb2TOF.dat
%     xdata_Solid = 1:dTOF1:1000; %[ns]
%     ydata_Solid = zeros(size(xdata_Solid));
%     ydata_Solid1 = dlmread(['KRb2TOF.dat'],'\t');
%     ydata_Solid1 = ydata_Solid1/max(ydata_Solid1)*ymaxKRb2;
%     xdata_Solid = xdata_Solid + (xdataKRb2-dTOF/2)*1000; %[ns]
%     ydata_Solid(xdata_Solid>36325&xdata_Solid<=36325+length(ydata_Solid1)*dTOF1) = ydata_Solid1;
%     plot(xdata_Solid/1000, ydata_Solid,'-r','LineWidth', 2.0);
%     hold off  
    
end
figure(h5);
% savefig(h, ['plots/summaryplot_K2Rb2_spectroscopy.fig']);
% saveas(h7, ['plots/TOFlineshape.png']);
% saveas(h3, ['plots/VMI_285nm_All.png']);